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INTRODUCTION 

In  this  report,  the  problem  of  estimating  the  three  parameters  of  the 
Carleton  damage  function  is  discussed.  The  technique  of  regression  analysis 
is  used  to  fit  the  damage  function  to  a  probability  of  kill  (&k)  matrix  with 
the  use  of  the  Statistical  Analysis  System  (SAS)  language  [3l,r4l  ,[5l.  The 
principal  output  of  the  study  are  the  estimates  Dg,  Rx ,  and  in  the  damage 
function  as  well  as  their  variance  -  covariance  matrix,  computed  for  a  set  of 
weapon/target  situations. 

The  converse  problem  is  also  addressed.  For  ?  target  located  at  (x,y) 
the  estimated  as  given  by  the  Carleton  damage  function  is  computed  as  well 
as  a  measure  of  the  error  in  Pk. 

The  report  is  divided  into  11  sections.  Section  II  gives  the  objective 
of  the  study.  Section  III  describes  the  matrix  while  Section  IV  digresses 
on  tne  three-parameter  Carleton  damage  function.  In  Section  V  and  VI  the 
various  methodologies  for  estimating  the  weapons  parameters  are  discussed. 
These  include  the  Roger's  (Snow)  ellipse  method,  the  linear  regression  method, 
and  the  non-linear  regression  method.  Sections  VII  and  VIII  are  related  to 
the  computational  phases  of  the  study  (SAS  language,  programs,  subroutines, 
etc.).  The  numerical  results  are  exhibited  in  Sectinn  IX.  In  Section  X,  the 
Taylor's  series  estimation  method  is  used  to  estimate  P^  and  to  compute  the 
error  on  p^.  Finally  Section  XI  is  a  concluding  section  with  recommendations. 
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SECTION  II 


OBJECTIVE  OF  THE  STUDY 

The  purpose  of  this  Study  is  to  determine  estimates  of  the  parameters  of 
the  generalized  three-parameter  Carleton  damage  function  as  well  as  the  error 
bounds  on  these  estimates.  The  mathematical  representation  of  this  function  is 

Pk  =  °0  expt"M^2  + 

x  y 

It  is  assumed  that  this  function  can  approximate  the  value  of  the  probability 
of  kill  appearing  as  elements  of  a  matrix  obtained  experimentally  as  a  result 
of  weapon  fragmentation  analysis  and  target  vulnerability. 

A  non-linear  regression  scheme  is  used  to  provide  a  best  fit  for  the 
Carleton  damage  function.  The  output  of  the  Study  consists  in  the  following: 

1.  The  estimates  of  the  three  parameters  Dq,  Rx,  and  Ry. 

2.  The  variance-covariance  matrix  of  Dq,  Rx,  and  Ry. 

The  estimates  obtained  agree  closely  with  the  estimates  using  Roger's 
Ellipse  Program.  Since  the  mathematical  basis  of  Roger's  program  cannot  be 
identified  and  since  the  program  does  not  generate  the  statistical  errors  on 
the  estimates,  it  was  necessary  to  use  a  non-linear  regression  technique  to 
arrive  at  these  errors. 

Using  the  Carleton  damage  function  and  its  parameter  estimates,  it  is 
possible  to  determine  for  any  weapon/target  system  the  probability  of  kill. 

Pi;,  at  any  location  (x,y)  given  that  the  weapon  bursts  at  (0,0).  In  addition, 
using  the  methodology  developed  in  [6],  the  variance-covariance  matrix  of 
Rx,  and  Ry  can  he  used  to  obtain  the  error  on  the  Pj,  value. 
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Parameter  estimations  were  performed  on  a  total  of  13  weapon/target  data 
sets.  Some  of  the  data  sets  were  obtained  using  the  same  veapon/target 
system,  but  using  different  weapon  impact  angles.  Different  impact  angles 
produce  different  P^  matrices,  thus  resulting  in  a  Carleton  damage  function 
with  different  parameters. 
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SECTION  III 
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DESCRIPTION  OF  THE  Ph  MATRIX 


A  ■ 

.  * 

■ 


t 
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The  matrix  is  an  array  of  mmbers  such  that  each  element  of  the  matrix 
provides  the  probability  of  kill  P^  for  a  given  point  target  positioned  at  a 
given  location  relative  to  the  point  at  which  the  weapon  bursts.  It  is 
assumed  that: 

1.  the  weapon  bursts  at  (0,0); 

?.  the  direction  of  the  x-axis  is  the  direction  of  range  (vertical 
direction) ; 

3.  the  direction  of  the  y-axis  is  the  direction  of  deflection  (horizontal 
direction) . 

The  maximum  size  of  the  Pk  matrix  printout  is  40  x  40  with  a  total  of  1600 
data  entries.  All  the  entries  are  non-negative  with  values  less  than  or  equal 
to  one.  Some  of  the  entries  are  simply  zero.  The  columns  of  the  array  are 
headed  by  the  ordinate  distance  between  the  point  of  weapon  burst  and  the 
horizontal  position  of  the  center  of  the  target.  These  distances  are  measured 
in  feet,  and  the  difference  in  distances  between  successive  columns  is 
constant.  Similarly,  the  rows  of  the  array  are  headed  by  the  abscissa 
distance  between  the  point  of  weapon  burst  and  the  vertical  position  of  the 
center  of  the  target.  Again,  these  distances  are  measured  in  feet,  and  the 
difference  in  distances  between  successive  rows  is  constant. 

As  would  be  expected,  target  points  which  are  close  to  the  weapon  burst 
center  have  a  high  probability  of  kill.  As  the  target  points  recede  from  the 
center,  the  probability  of  kill  decreases  until  it  reaches  zero.  It  is 
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conceivable  that  even  at  the  point  £0,0)  the  probability  of  (till  could  be  less 
than  one;  this  is  usually  a  characteristic  of  the  target/weapon  system  un^er 
consideration. 

In  general,  the  Pj<  matrix  will  combine  the  effect  of  weapon  fragmentation 
ano  blast.  For  the  study  under  consideration,  only  fragmentation  effects  are 
accounted  for. 

On  an  actual  matrix,  only  the  upper  right  and  the  lower  right  entries 
appear.  The  upper  left  and  the  lower  left  entries  do  not  appear  as  it  is 
usually  assumed  that  the  effect  of  the  weapon  burst  is  symmetrical  anout  *;ne 
range  or  x-axis. 

In  general,  if  the  weapon  impacts  at  90  degrees,  the  Pk  matrix  is 
symmetrical  about  the  deflection  or  y-axis.  If,  on  the  other  hand,  the  weapon 
impacts  at  an  angle  less  than  90  degrees,  considerable  asymmetry  may  result 
between  the  upper  right  entries  and  the  lower  right  entries  of  the  Pj<  matrix. 
The  full  display  of  an  equiprobabil ity  curve  as  appearing  in  the  P^  matrix 
resembles  the  wing  contour  of  half  a  butterfly. 
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SECTION  IV 


THE  CARLETON  DAMAGE  FUNCTION 

The  Carleton  damage  function  is  a  mathematical  formula  which  gives  the 
probability  of  kill  P|<  of  a  weapon  bursting  at  the  point,  (0,0)  upon  a  target 
located  at  (x,y).  Both  weapon  and  target  are  idealized  as  points,  and  the  only 
prevailing  effect  is  due  to  weapon  fragmentation.  The  formula  is: 

Pk  *  0o  **Pl-0o(~*  ^)j 
x  y 

In  its  more  general  form,  the  function  has  three  parameters: 

1.  The  parameter  Rx  which  is  the  weapon  radius  in  the  direction 
of  range  or  x-axis  and  is  measured  in  feet. 

?.  The  parameter  Ry  which  is  the  weapon  radius  in  the  direction 
of  deflection  or  y-axis  and  is  measured  in  feet. 

3.  The  parameter  D0  *hich  represents  the  maximum  probability  of  kill 
which  occurs  at 

Of  course,  nn  must  always  be  a  positive  quantity  whose  value  is  less  than  or 
n-iml  to  one.  For  most  of  the  weapons  under  consideration,  Ry  >  RK ,  since  the 
fragmentation  effect  in  the  direction  of  deflection  is  significantly  more 
marked  than  in  the  direction  of  range. 

In  general,  the  parameters  of  the  Carleton  damage  function  will  depend 
on : 

1.  The  impact  angle  of  the  weapon  (in  the  direction  of  range). 

?.  The  weapon  characteristics. 

3.  The  target  characteristics. 

A.  The  weapon  velocity. 
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When  P^  is  plotted  against  x  and  y,  the  shape  of  the  generated  surface  is 
that  v.f  a  bivariate  Gaussian  distribution.  The  equiprobabi 1 ity  contour  lines 
are  concentric  ellipses  with,  in  general,  the  major  axis  along  the  y-axis  and 
the  minor  axis  along  the  x-axis. 

The  theoretical  values  of  obtained  from  the  damage  function  are  such 

that  they  are  always  positive  for  any  (x,y)  values.  That  is,  it  is  assumed 

that  Pt<  is  always  positive  irrespective  of  the  magnitudes  of  x  and  y.  In 

practice,  however,  the  coverage  area  of  the  weapon  due  to  fragmentation  is 

finite  and  outside  of  a  certain  contour  about  the  origin,  P^O. 

The  theoretical  equi probability  contour  lines,  which  are  concentric 

ellipses,  approximate  the  experimentally  obtained  contour  lines  which  are,  in 

general,  butterfly  shaped.  Thus,  in  using  the  obtained  from  the  Carleton 

function  at  a  given  (x,y),  a  discrepancy  will,  in  general,  result  between  the 

true  mean  value  of  P^  (the  observed  sample  mean)  and  the  calculated  mean  value 

Gf  Pk  which  is  the  theoretical  E[Pk]  from  the  model.  This  difference  is 

reflected  in  terms  of  a  standard  deviation  ap  .  Thus,  the  expression 

k 

ECPJ  t 

K 

will  state  that  in  at  least  75  percent  of  the  cases  (using  Chebyshev's 

theorem)  the  true  mean  value  of  Pj.  will  be  within  the  above  range,  a  can  be 

k 

calculated  (see  Section  X)  explicitly  as  a  mathematical  formula  involving  the 

estimates  of  Oq,  Rx,  and  Ry  as  well  as  their  variance-covariance  matrix.  The 

numerical  value  of  o  can  then  be  obtained  using  the  outputs  of  the  non- 

k 

linear  regression  analysis  which  provides  the  numerical  values  of  Dq,  Rx,  and 
R y  as  well  as  their  variance-covariance  matrix. 
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SECTION  V 


ESTIMATION  OF  THE  WEAPONS  PARAMETERS  (1) 

In  this  section  Roger's  ellipse  program  Is  briefly  discussed  and  a 
summary  of  the  linear  regression  procedure  is  presented.  The  non-linear 
regression  procedure  based  on  the  Gauss-Newton  method  is  discussed  in  detail 
in  Section  VI. 

1 .  Roger's  Ellipse  Program 

To  estimate  the  parameters  Dq,  Rx  and  Ry  of  the  Carleton  damage  function, 
Mr.  Roger  Snow,  formerly  with  the  RANO  Corporation,  developed  an  algorithm 
which  appears  in  a  program  listing  as  a  subroutine  under  the  heading  "Roger's 
Ellipse  Program,"  Unfortunately,  the  subroutine  is  undocumented,  and  it  was 
not  possible  to  derive  the  methodology,  the  mathematical  rationale,  and  the 
equations  leading  to  the  algorithm.  The  program  subroutine  provides  a  single 
estimate  number  for  Dq,  Rx,  and  Ry.  Despite  the  fact  that  the  Carleton  damage 
function  specifies  a  value  of  Dq  in  the  range  0  <  Oq  £  1,  the  algorithm  has 
sometimes  generated  a  value  of  Dq  >  1,  which,  in  effect,  means  that  the 
largest  probability  of  kill  obtained  analytically  may  take  on  values  greater 
than  one.  Finally,  the  ellipse  program  does  not  compute  the  errors  associated 
with  the  estimation  of  the  three  parameters.  (See  Appendix  A  for  a  listing  of 
the  program.) 
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2.  Linear  Regression  Method 

Recall  that  the  expression  for  the  Carleton  damage  function  is 

pk » o0  mp(-o0c4  *  4] 

Rx  y 

By  taking  natural  logarithms  on  both  sides  of  this  expression  one  obtains 

2  ^02 

In  Pfc  =  In  Oq - ^  x  -  -j  y 

R  R 

x  y 

This  can  be  made  a  linear  expression  in  u=x£  and  v=y  and,  hence.  In  Dq, 

2  2 

(-  D0/Rx)  and  (-  DQ/R^}  can  be  estimated  using  the  classical  techniques  of 
linear  regression.  Details  of  the  methodology  when  Dq=1  is  given  in  Tfil. 

Unfortunately,  the  technique  is  ill-suited  to  provide  accurate  estimates 
of  the  parameters  under  consideration.  There  are  three  main  factors  which 
contribute  to  this  situation; 

a.  Since  Pk  lies  in  the  range  0  <  Pk  <  1,  In  Pk  becomes  very  sensitive  to 
small  variations  in  Pk. 

b.  Data  entries  for  which  Pk=0  In  the  matrix  must  be  discarded  since 
they  generate  infinite  values  for  In  Pk. 

c.  Points  which  are  far  from  the  origin  are  weighted  more  strongly  than 
are  those  points  near  the  origin.  Hence,  a  poor  fit  near  the  origin 
results. 

To  see  this,  note  that  the  sum  of  squares  of  the  residuals  for  the 
linearized  model  has  the  form 

l  0"  Pk  -  In  P*  f 
i=l  Ki  i 
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where  n  =  number  of  observations; 

P.  *  calculated  value  of  the  1^  data  point; 

P£  =  experimental  measured  value  of  the  I**1  data  point. 


m 


then  note  that 


6  Pk  “  Pk  -  Pk 
i  i  K1 


1n  pk  "  ln  pk  =  1n  (“*4 

i  K1  Pk 
*1 


»  ln(  - 


V  +  6  Pk: 


=  ^  (1+  tt~  •) 

P.. 


If  the  6  P,  ‘s  are  small,  then 
Ki 


&  P.  6  P. 

k,  k . 


Hence,  the  sum  of  squares  of  tre  residuals  is  given  approximately  by 


»  k 

l  (-r-42 

i=l  P, 


On  the  other  hand,  the  sum  of  squares  of  the  residuals  for  the  non-linear 
regression  model  is  simply 


i  («  »k  r 

1=1  *i 
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*  2 

The  (l/P,  )  factor  in  the  linearized  model  acts  as  a  weighting  factor  which 

i 

gives  additional  weight  to  those  points  for  which  P.  is  small,  i.e.,  far  from 
the  origin. 

The  main  advantage  of  the  method  is  that  analytic  techniques  are 
available  to  obtain  closed  form  expressions  for  the  estimates,  as  well  as 
their  variance-covariance  matrix.  A  more  accurate  technique  is  that  of  non¬ 
linear  regression  which  will  be  discussed  later. 

The  linear  regression  method,  however,  may  be  used  to  obtain  an  initial 
estimate  of  the  parameters  which  then  can  be  used  as  a  starting  point  for  the 
non-linear  regression  technique.  To  simplify  the  computational  routine,  an 


SECTION  VI 


ESTIMATION  OF  THE  WEAPONS  PARAMETERS  (2) 

In  this  section,  a  non-linear  regression  procedure  based  on  the  Gauss- 
Newton  technique  is  discussed  (see  e.g.,  [1],[?D. 

1 .  Non-Linear  Regression  Method 

Suppose  that  it  is  required  to  fit  a  function  of  the  form 

Pi  =  M(xi  ,y i;  A,B,C)  (1) 

i=l,2,...,n,  to  a  large  collection  of  n  data.  Here  x,  and  yj  are  a  pair  of 
independent  variables,  and  is  the  corresponding  dependent  variable.  A,B, 
and  C  are  the  parameters  of  the  model  whose  values  are  to  be  estimated  through 
a  least  square  procedure. 

The  residual  at  each  point  i,  i=l,?,...,n,  is  defined  by 

ri  “  M(x.,yi;  A,B,C)  -  Pi  (2) 

Tnen,  tm?  parameters  A,B,  and  C  are  determined  by  minimizing  the  sum  of 
squares  of  the  residuals.  For  mathematical  convenience,  the  objective 
function  adopted  has  the  form 

f(A,B,C)  =  i  [  r?  (3) 

i  =  l 

To  minimize  this  function,  set  its  partial  derivatives  with  respect  to  A,R, 
and  C  to  zero  and  thus  obtain 


(4) 


r,  3r. 


n  3r  = 


n  ar. 

l  riir-n 

i-i 


Using  (2)  and  for  brevity,  setting  M,  equal  to  M(  x-j  ,y^  ,A,R,C) ,  equations  (A), 


(5),  and  (6)  become 


n  3M . 

I  (Mi  -  Pi  I  5T  =  n 

i=l 


n  3M, 

I  ("i  -Pi)ar=° 

1  =  1 


n  3M. 

I  (Mi  -  P,)  ar"  =  o 
1*1 

These  last  three  equations  must  now  be  solved  for  A,B,  and  C,  Generally, 
these  equations  non-linear  and  a  standard  method  of  solution  is  the 

iterative  method  of  Newton. 

For  the  moment  define 


F(A,B.C)  =  I  (Mi  -  p,0  jjs^- 
1=1 


II 

G(A.B.C)  =  I  (Mi  -  Pi)  sg- 

i  =  l 
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(12) 


n  3M. 

H(A,B,C)  *  \  (M.  -  p^)  yjt— 

Let  A*,  8*,  and  C*  be  the  solution  to  the  simultaneous  set  of  equations 


F(A,B,C)  =  0 

(13) 

G(A,B,C)  =  0 

(14) 

H(A,B,C)  =  0 

(15) 

Also, let  A0,  B°, and  C°  be  an  approximate  solution  of  (13)  and  (14)  and  (15), 
while  A4" ,  B+,  and  C+  represent  an  approximate  solution  to  (13),  (14),  and, 
(15)  which  is  an  improvement  over  A0,  8°,  C°. 

By  Taylor's  series 

F(A*.  B*.  C*)  *  0  =  F(A°,  8°,  C°)  +  (A*  -  A0) 


+  <B*  -  B°)  ff  +  (C  -  C°)  ... 


,0,  3F 


(16) 


G( A* ,  B*.  C*)  =  0  =  G( A°,  Bu,  Cw)  +  (A*  -  A0)  || 


+  (B*  -  B°)  ||  +  (C* 


C°)|§  + 


(17) 


H(A*.  e*.  C*)  -  0  =  H(A°,  B°,  C°)  +  (A*  -  A°) 


+  (B*  -  B°)  ||+  (C* 


(18) 
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where  all  the  partial  derivatives  are  computed  at  the  point  (A0,  B°,  C°) .  If 
the  right  sides  of  (16),  (17),  and  (18)  are  retained  as  infinite  series,  one 
could  solve  in  principle  for  A*.  B*,  and  C*.  However,  if  the  series  are 
truncated  just  before  the  second  derivative  terms,  the  equations  cannot  be 
solved  for  A*,  8*,  and  C*,  but  instead  may  be  solved  for  A+,  B+,  and  C+.  This 
is  easily  done  since  the  resulting  equations  are  linear.  Writing  the 
truncated  versions  of  (16),  (17),  and  (18)  in  matrix  form,  one  obtains  at  the 
point  (A°,B°,C°) 


3F 

3F 

3F 

3  A 

38 

3C 

3G 

3G 

3G 

3  A 

38 

3C 

3H 

3H 

3H 

3  A 

3B 

3C 

A+  -  A0 


B+  -  B° 


C+  -  C° 


L 


By  means  of  (10),  (11),  and  (12),  the  quantities  F,  G,  and  H  may  be  removed; 
consequently,  the  matrix  on  the  left-hand  side  of  (19)  may  be  written  as  the 
sum  of  two  matrices,  namely. 


3F 

3F 

3F 

u 

F5 

37 

3G 

3G 

3G 

3A 

3B 

3C 

3H 

3H 

3H 

TK 

3& 

37 

111  ™  =  s  ♦  T 
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where 


1  =  "  (M1  -  Pi) 

i=l 


and 


32M4 

a2^ 

3^ 

lev; 

< 

1 

ala" B 

3A3C 

3 

A 

3B3A 

3B? 

3B3C 

A 

aV 

1 

3^ 

3T!3A 

3C3B 

ac  = 

3Mi 

3Mi 

3Mi 

3M. 

1 

3Mi 

3M. 

TT 

aT' 

alT 

38“ 

al~ 

ac 

3M. 

l 

3M- 

1 

3M. 

3M. 

3Mi 

aM 

i 

a¥ " 

3A  * 

38" 

aT 

36 

a(T 

3Mi 

3Mi 

3Mi 

3Mf 

3Mi 

3Mi 

aTT 

3A" 

alT 

3B“ 

ac 

an- 

The  ngnt-hand  side  of  expression  (19)  may  be  written  as 


F 

3"i 
a  A 

n 

3M. 

6 

=  1  -  Pi ) 
i  =1 

‘ \ 

TB 

H 

3Hi 

“ac 

. 

. 

In  terras  of  S,  T  and  b,  (19)  becomes 


(21) 


(22) 


(23) 
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A*  -  A0 


(s  +  D 


8+  -  B° 
C+  -  C° 


=  -  b 


After  solution  one  gains 


“ 

A+ 

A0 

B+ 

If 

B° 

C+ 

c°  ! 

i 

- 

L  J 

(S  +  I)'1  b 


(24) 


(25) 


where  it  is  understood  that  the  right-hand  side  of  (25)  is  evaluated  at  the 
point  (A°,B°,C°). 

From  (25),  one  obtains  A+,  B+,and  C+  which  are  now  redefined  as  A0,  B°, 
C°,  and  the  process  is  repeated  until  the  vectors  [A0,  B°,  C°]  and  [A+,  B+, 
C+]  are  the  same,  or  at  least  until  they  become  arbitrarily  close  to  each 
other.  The  amount  of  work  is  staggering,  and  simplifications  are  sought. 


2.  Case  when  M(x,y;  A,B,C)  is  linear 

If  M(x,y;  A,B,C)  is  linear  in  A,B, and  C,  then  two  important 
simplifications  result. 

a.  becomes  zero,  and 

b.  Mi  i  Mt  (Xi,yi;  A°,B°,C°) 

3M. 


3M. 

-  A°  3A1  +  VT  +  c  3C 


„  ....  SM . 

0  1  „0  1 


(26) 


Setting  to  zero  and  multiplying  by  ^transforms  (25)  into: 

17 


.  *~irTT  1 1  in 


rftjM  ii  i  *’»»»  ■  aAjij 


For  the  linear  case,  because  of  (26),  this  is  the  same  as 


r  ] 

A° 

5i r 

n 

3M. 

B° 

=  l  M1 
i  =  l 

l 

3B~ 

3M- 

C° 

1 

J 

L  ^  . 

(30) 


Hence,  the  first 


two  terms 


on  the  right 


side  of  (27)  cancel  yielding 


ri 

a, p< 


3M. 

1 

aTT 


3M. 

1 

3B~ 


5£~ 


1 


(31) 


Thus,  A+,  B+,  and  C+  may  be  solved  for  without  iteration.  These  values  are  in 
fact  A*,  B*  and  C*. 


3.  The  Gauss- Newton  Method 

Returning  to  the  case  when  M(x,y;  A,B,C)  is  not  linear  in  A.  B, and  C,  one 
may  use  (25),  in  the  iterative  fashion  outlined,  in  order  to  approximate  A*, 
B*,  and  C*.  This  method  is  tedious  and  has  been  called  the  full-Newton  method. 
As  the  full-Newton  method  proceeds,  b  becomes  closer  and  closer  to  zero,  from 
(21),  one  may  see  that  is  linear  in  the  residuals  (M^-p^).  Since  it  is 
reasonable  to  suppose  that  these  •■esiduals  decrease  during  the  iteration 
process,  it  follows  th3t  becorMS  less  and  less  important  as  the  iterations 
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evolve,  The  Gauss-Newton  method  exploits  this  notion  by  setting  equal  to  ’  •  -j 

zero.  Hence,  the  Gauss-Newton  method  uses  the  formula  J 


- 

r 

-  y  v-j 

A+ 

A° 

-  r1  b 

B+ 

(32) 

-  •  i 

C+ 

C° 

■J 

- 

J 

in  an  iterative 
right  hand  side 
C°  are  replaced 
continues  until 
to  each  other, 
derivatives. 


fashion.  That  is.  A0,  B°  and  C°  are  used  to  evaluate  the 
of  (32)  and,  hence,  to  calculate  A+.  B+  and  C+.  Next,  A0.  B°  and 
by  the  A+,  B+  and  C+  that  were  just  computed.  The  process 
the  vectors  [A+,B+,C+]  and  [A°,B°,C0]  become  arbitrarily  close 
The  Gauss-Newton  method  does  not  involve  the  use  of  second 


4.  itenarr 

Before  passing  on  to  an  example,  let  it  be  noted  that  some  authorities 
define  the  matrix 

3M1  aMj 

I  3A~  3B~  3C  I 


3Mn 
3  A 


3Mr 

3B 


3«n 

3C 


(33) 


From  [??).  X  is  easily  computed  as 


r  =  ,T  , 

•  V  W 


(34) 
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Also,  if  one  defines  the  vectors 


one  immediately  obtains  from  (23) 

b_  =  JT  f.M  -  £l 


(36) 


(36) 


The  vector  [M_ -  jj]  is  called  the  residual  vector.  The  root  mean  square 
( RMS)  associated  with  the  residual  vector  is  defined  by 


where  f  is  the  value  of  the  objective  function  given  in  (3)  and  n  is  the 
number  of  observation  in  the  data  set.  The  RMS  is  an  aggregate  measure  of  the 
residuals,  and  as  the  iterations  evolve,  the  value  of  the  RMS  should  decrease 
and  converge  to  a  fixed  value. 


5.  Example 


Let  the  function  in  question  which  is  to  be  used  to  fit  the  data  be 


pk  =  °o  expt'Mr? + 

Rx  Ry 


Let  the  given  data  be 


(38) 


i 

xi 

yi 

\ 

i 

0 

0 

.9 

2 

0 

1 

.5 

3 

1 

0 

.6 

4 

1 

1 

.3 

The  parameters  to  be  estimated  are  Og,  Rx  and  Ry.  To  obtain  the  initial 
values  of  the  parameters,  namely  Dg,  R°  and  R°,  use  a  linear  model.  To  this 
end,  take  the  natural  logarithm  of  (38)  and  obtain: 

Cg  2  ^0  2 

In  Pk  -  In  Og  -  -g  x2  -  J  yZ  (39) 

Kx  Ry 

An  initial  value  of  Dg  may  be  taken  as  Dg=l.  However,  this  is  not  necessary 
as  will  be  noted  in  the  sequel.  Let 


a  =  In  Dg, 


(40) 
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From  (35) 


■ 

-» 

‘ 

Pi 

1  n 

Pk 

K1 

pk 

-.1053605 

P2 

In 

-.6931472 

= 

2 

= 

P3 

In 

pk. 

-.5108256 

3 

P4 

In 

\ 

-1.2039728 

. 

. 

.  • 

The  right-hand  side  of  (31)  is 


-2.513306124 

-1.714798428 

-1.897119985 


Hence,  (31)  becomes 


- 

■ 

“  ■« 

4  2  2 

+ 

a 

-2.513306124 

2  2  1 

a+ 

-1.714798428 

2  1  2 

b+ 

-1.897119985 

- 

. 

The  solution  is 


a+  =  -.0790203868 
a4-  =  -.4581453659 
b+  =  -.6404669227 
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From  (41),  (42),  and  (43)  it  follows  that 


Dq  =  .9240210864 
R°  =  1.420166592 
R°  =  1.201137119 

We  first  compute  the  RMS  value  for  the  non-linear  function  (38).  The 
is  obtained  by  substituting  (45)  in  (38)  and  calculating  it  for  each  i 
to  yield 

- 

.9240210864 

.4870018732 

M  = 

.5844022479 

.3080070289 

m 

The  p_  vector  is  immediately  given  from  the  data  set.  The  logarithmic 
transformation  is  not  used  for  the  non-linear  function. 


(45) 


vector 
1,2, 3, 4 
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The  residual  vector  for  (D~,  R°  R°)  is 

'J  x  y 


.0240210854 


M  -  £ 


-.0125981268 

-.0155977521 


.0080070289 


This  results  in  a  RMS  value  of  .0162278023.  Values  in  (45)  will  now  be  used 
as  initial  values  in  the  Gauss-Newton  method  where  the  non-linear  function 
(38)  is  being  employed.  Note  that 


3M  3Pk  .  Y2  v2  Y2  v2 

arST  =  3d7  =  t1  '  °o^r?  +  „2^  expt‘  °o^72  +  nx^l 


RX.  Ry 


Ry 


3M  3Pk  9  Dn  v2  J 

3fT=3R-=2x  3“expf-  Dobf  + 

X  x  Rx  Rx  Ry 


aM  3pk  ?  Dn  y2  y2 

4-e*P[-Dn(V^)] 


3R  3R 

y  y 


0  R2  R2' 

Rx  Ry 


The  new  J_ matrix  as  defined  in  (33)  is  now  formed  using  (461 ,  (47),  and  (48), 
all  evaluated  a.  the  point  (D°,  R°,  R°)  as  given  by  (45). 


j  .1894905697  0 

I  .3426989609  .3 


-.0328707630  1987259717 


.5193555108 


.3770560205  0 


.3284692657 
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Using  (52)  a  new  and  £  vector  can  be  calculated  resulting  in  a  residual 
vector  of 

.0060475147 

„  -.0117137842 

M  -  £  = 

-.0102837618 

.0178059923 

The  corresponding  RMS  is  .0122138516. 

Thus,  one  iteration  of  the  Gauss-Newton  method  has  been  completed.  As 
can  be  seen,  the  RMS  value  associated  with  the  residual  vector  decreased  from 
.0162278023  to  .0122138516  which  indicates  that  the  method  is  converging. 

6.  The  Variance  -  Covariance  Matrix 

Whereas,  in  the  linear  least-square  model  the  variance-covariance  matrix 
of  the  parameter  estimates  was  uniquely  defined  and  could  be  derived 
analytically,  the  situation  is  quite  different  for  non-linear  least  squares. 

In  non-linear  least-squares  several  corresponding  variance-covariance  matrices 
are  used,  all  of  them  being  approximations.  In  addition,  they  are  not  as 
reliable  as  the  linear  least  squares  one.  !n  fact,  some  authors  (see  e.g., 

Cll)  even  caution  users  of  non-linear  least-squares  software  about  an 
indiscriminate  use  of  the  variance-covariance  estimates. 

n 

One  particular  variance-covariance  matrix  that  is  often  used  is  (RMS)  _T_f , 
where  J[j  is  the  final  value  of  the  matrix  (34)  computed  at  the  final  estimate 
values  of  the  parameters,  and  (RMS)  is  the  root  mean  square  associated  with 
the  residuals. 
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SECTION  VII 


COMPUTATIONAL  ASPECTS  OF  THE  PROBLEM 

1-  Program  Development 

During  the  development  of  the  programs,  factors  that  were  considered 

were: 

( a)  Accuracy  of  results 

Would  significant  error  result  due  to  roundoff,  particularly  due 
to  the  large  number  of  complex  mathematical  operations  required  on  a  large 
data  set?  Conceivably,  a  small  roundoff  error  occurring  at  the  beginning 
could  propagate  and  result  in  major  errors  in  the  estimates  of  the  parameters. 

(b)  Simplicity  and  flexibility 

The  procedures  developed  should  be  easy  to  understand  and  be  adaptable 
to  several  computer  facilities. 

(c)  Data  management  and  additional  analysis  features 

Adequacy  of  data  storage  should  be  a  prime  requirement.  Also,  the  pro¬ 
grammer/analyst  might  need  to  perform  additional  analyses,  plots,  charts,  etc. 

With  these  considerations  in  mind,  it  was  decided  to  use  the  Statistical 
Analysis  System  (SAS)  language  to  perform  the  desired  linear  and  non-linear 
regression  analyses.  For  detailed  description  the  reader  is  referred  to  [5l. 

For  the  linear  regression  case,  the  PROC  REG  procedure  was  used.  This 
procedure  is  based  on  the  least  squares  principle  to  produce  Best  Linear 
Unbiased  Estimates  (BLUE)  under  classical  statistical  assumptions.  For  the 
non-linear  regression  case,  the  PROC  NLIN  procedure  was  used.  Several  non¬ 
linear  estimation  methods  are  available  using  this  procedure  (e.g.,  Gauss- 
Newton,  Marquardt  and  steepest  descent). 
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2.  Program  Performance 


Two  sets  of  programs,  one  for  the  linear  case  and  the  other  for  the  non¬ 
linear  case,  were  developed.  To  check  the  degree  of  roundoff  error,  a  15- 
observation  data  set  was  analyzed  using  the  PROC  REG  procedure.  The  parameter 
estimates  and  computed  variances/covariances  were  in  agreement  with  hand- 
calculated  results  to  the  eighth  decimal  place.  The  programs  were  then  tested 
on  a  full  size  { 1 ,600-observation)  data  set  and  produced  output  in  the  desired 
format . 

Two  very  short  subroutines  were  written  and  used  to  plot  the  residuals 
(observed  P^  -  predicted  Pk)  versus  the  x(range)  values  or  the  y  (deflection) 
values  and  to  provide  a  printout  of  the  residuals.  The  plotting  routine  was 
used  extensively  in  evaluating  output  result. 

For  the  PROC  NLIN  procedure,  both  the  Gauss-Newton  and  the  Marquardt 
methods  were  used  on  a  full-size  data  set  with  results  in  agreement  to  the 
third  decimal  place  for  the  parameter  estimates  and  the  eighth  decimal  place 
for  the  variance  -  covariance  matrix.  The  Gauss-Newton  method,  as  described 
in  Section  VI,  was  the  one  that  was  arbitrarily  selected  and  used  on  the  rest 
of  the  data  set. 

Finally,  a  number  of  runs  under  the  PROC  NLIN  procedure,  using  arbitrary 
extreme  starting  values  for  the  parameters  (PARKS  statement),  indicated  that 
close  starting  values  for  the  Rx  and  Ry  parameters  were  not  necessary  for 
accurate  results,  but  might  lead  to  significant  savings  in  computer  time. 
However,  different  starting  values  for  the  Dg  parameter  sometimes  resulted  in 
a  non-optimal  set  of  parameter  estimates.  Hence,  the  PARKS  statement  in  the 
PROC  NLIN  procedure  is  written  to  assure  a  close  starting  value  for  Dg  in  the 
iterative  phase. 
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A  question  addressed  was  the  propriety  of  using  the  PROC  REG  procedure 
under  the  assumption  that  the  two-parameter  Carleton  damage  function  (Qq*1) 
was  a  reasonable  model,  and  that  the  function,  as  previously  discussed  in 
Section  v,  can  be  linearized  in  order  to  estimate  Rx  and  Ry.  To  answer  this 
question,  the  procedure  was  used  on  a  full  size  (l,60Q-observ3tions)  data  set 
as  provided  by  one  of  the  matrices.  A  printout  and  plot  of  the  residuals 
obtained  showed  that  the  residuals  ranged  between  ±0.?  and  were  particularly 
large  (in  absolute  value)  in  the  area  of  weapon  impact  (origin). 

Using  the  same  data  set,  PROC  NUN  procedure  was  used  under  the 
assumption  that  the  three-parameter  Carleton  damage  function  (0  <  DQ  <  1)  was 
a  better  model.  In  this  instance,  there  was  no  need  to  manipulate  the  model 
as  in  the  previous  case  or  to  discard  values  of  Pk=0.  The  residuals  obtained 
ranged  between  ±0.15  in  the  area  of  the  origin  and  were  much  better  elsewhere. 

The  poor  fit  in  the  neighborhood  of  the  origin  is  due  to  the  preponder¬ 
ance  of  the  data  points  being  outside  that  neighborhood.  The  points  in  the 
neighborhood  of  the  origin  are  simply  outnumbered  by  the  other  points.  The 
values  of  near  the  origin  have  almost  no  influence  on  the  values  of  the 
parameters  of  the  model . 

For  comparison  purpose,  a  linearized  regression  analysis  was  performed  on 
the  data  set  corresponding  to  the  following: 

Weapon:  5EAKL,  Target  8T1.5.  Impact  angle  15°  Velocity  900  ft/sec. 

Using  the  three-parameter  Carleton  damage  function.  The  residual  sum  of 
squares  (RSS)  obtained  for  the  P^  values  was  14.P0Q.  This  was  considerably 
more  than  the  RSS  obtained  using  the  non-linear  regression  mod<»l ,  which  was 
The  estimates  of  0g,  Rx,  and  Rv  for  the  two  models  were  as  follows: 


3 -parameter 

1  inear  model 

3-parameter 
non-linear  model 

"o 

0.020221 

0.75515 

Rx 

99.570 

72.9501 

Ry 

104.332 

269.2381 

RSS 

14.209 

3.269 

In  the  least-square  sense,  the  non-linear  regression  model  yields  much 

tetter  parameter  estimates  than  the  linearized  regression  models.  Also,  as 

previously  stated,  the  two-parameter  linearized  model  yielding  estimates  for 

R  and  R  obtained  from  PROC  REG  procedure  were  used  as  the  initial  values  for 

*  y 

PROC  NUN  procedure. 


! 


I 


I 
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SECTION  VIII 


COMPUTATIONAL  ROUTINE 


1.  Introduction 

In  this  section,  the  computational  routine  associated  with  SAS  is 
detailed. 

The  primary  input  that  is  provided  by  the  Air  Force  consists  of  the  40  x  40 
probability  of  kill  (P^)  matrix  for  a  given  weapon/target  combination  speci¬ 
fied  for  a  given  weapon  impact  angle  and  impact  velocity.  This  matrix  is 
given  as  a  computer  printout. 

2.  Oesign  Information 

A  copy  of  the  computer  program  used  appears  in  Appendix  P. 

In  general,  this  program  was  designed  to  output  all  the  desired  estimates 
using  both  the  linear  and  non-linear  regression  procedures,  with  the  user  only 
required  to  provide  the  input  data.  In  addition,  this  program  is  written  so 
that  the  input  data  set  needs  only  contain  3  simple  lines  of  code  and  the  1600 
probability  of  kill  (Pk)  values  typed  (or  punched)  ten  to  a  line  (or  card). 

The  majority  of  the  code  in  this  program  is  necessary  only  to: 

(1)  Manipulate,  generate,  label,  or  properly  format  the  input  and 
program  generated  data  sets  for  the  linear  and  non-linear  analyses. 

(2)  Generate  the  properly  formatted  output. 

Although  a  thorough  understanding  of  SAS  is  necessary  to  comprehend  these 
portions  of  the  code,  this  should  pose  no  problem  in  adapting  this  program  for 
execution  on  other  computers  or  at  other  facilities. 

It  should  be  noted,  however,  that  the  two  INFILE  statements  in  the 
program  are  used  to  read  the  input  data  set  from  a  disk  file  and  to  read  a 
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program  -  generated  external  disk  file.  Statements  must  be  included  in  the 
JCl  (for  a  batch  run)  or  on  the  terminal  (prior  to  an  interactive  run)  to 
associate  the  correct  disk  files  with  the  1NFILE  statements. 

3.  Input  Data 

The  input  data  set  may  be  entered  onto  a  computer  disk  file  or  onto 
punchcards  in  the  following  format,  each  line  representing  one  punchcard  or 
one  line  on  the  terminal: 


X 

19.8 

13. 

Y 

13.7 

0.0 

l 

0.0 

0.0 

m  # 

e 

.0011 

.0013 

.0017  .  .  . 

.  .  .  .0001 

.000? 

.0017 

.00?0  .  .  . 

.  .  .  .0007 

160  lines 


.  .0001  ,0007  .0008  . 0004 

10  values 

-  - » 

To  input  each  separate  data  set,  refer  to  the  example  above  and  do  the 

foil owi ng ; 

Line  a  1 

The  first  entry  is  always  the  letter  X.  The  second  entry  will  be  the 
difference  between  successive  range  values  in  the  data  set.  Note  that  this 
value  is  a  constant  for  each  data  set.  The  third  entry  will  be  computed  as 
fol lows: 
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most  negative  range  value _ 

difference  between  successive  1 

range  values 

which  is  always  an  integer. 

Line  #2 

The  first  entry  is  the  letter  Y.  The  second  entry  is  the  difference  between 
successive  deflection  values  in  the  dataset.  Note  that  this  value  is  a 
constant  for  each  data  set.  The  third  entry  will  be  the  number  0.0. 

Line  #3 

The  first  entry  is  the  letter  2 .  The  second  and  third  entries  are  both  the 
number  0.0 

Line  #4  -  Line  #163 

Refer  to  the  Pk  values  as  they  are  printed  on  the  Air  Force  computer  output. 
There  are  four  pages,  each  with  400  Pk  values.  Each  page  contains  the  Pk 
values  corresponding  to  10  different  deflection  values  and  all  40  range 
values.  Each  of  the  160  lines  corresponds  to  one  line  of  input  data,  going  in 
order  from  first  page  to  last,  top  to  bottom. 

4.  Program  Execution 

The  first  part  of  the  program  consists  of  the  code  necessary  to  read  in 
the  input  data  set,  generate  the  range  and  deflection  values  corresponding  to 
each  Pk  value,  and  transpose  all  these  values  into  the  proper  format  for  the 
analyses. 

The  second  part  of  the  program  starts  with  PROC  REG  statement  and  MODEL 
statement  necessary  to  call  the  PROC  REG  procedure  and  specify  the  model  to  be 
fitted,  respectively.  Following  these  are  several  more  statements  used  for 
manipulation  of  the  PROC  REG  output. 

35 


ft  i 


-w V- 

t.-.  ^ 


rr.  i|i 


The  next  section  begins  with  the  statement; 

PROC  PRINTTO  UNIT  =  20; 

(This  causes  the  PROC  NUN  non-linear  regression  output  to  be  routed  to  an 
external  disk  file  which  is  referenced  later  in  the  INFILE  FT20F001 
statement.)  Following  is  a  brief  description  of  the  statements  used  to 
perform  the  non-linear  regression: 

(1)  PROC  NUN  statement  followed  by  the  DATA  =  statement.  The  DATA  = 
statement  specifics  the  data  set  to  be  analyzed.  A  METHOD  *  statement  may 
also  be  included  specifying  the  mathematical  procedure  to  be  used.  The 
default  method  is  Gauss-Newton. 

(2)  PARMS  statement  followed  by  the  RX  =  ,  RY  =  ,  and  DO  =  assignments 
specifying  the  grid  of  initial  parameter  values  to  be  used.  The  procedure 
computes  the  sums  of  squares  for  error  for  each  possible  combination  and  uses 
the  parameter  values  in  the  grid  corresponding  to  the  least  such  sums  of 
squares  as  the  initial  parameter  values  in  the  iterative  procedure. 

(3)  ROUNDS  statement  followed  by  code  designating  the  regions  of 
consideration  for  each  parameter  estimate.  The  procedure  will  not  consider 
narumpter  values  outside  the  boundary  regions  during  the  iterative  process. 

(A)  MODEL  statement  designating  the  model  to  be  fitted,  followed  by 
statements  of  the  form  0ER.00  =  which  give  the  partial  derivatives  of  the 
model  expression  with  respect  to  the  parameters  to  be  estimated. 

Section  4,  beginning  with  the  statement 
DATA  TWO; 

consists  entirely  of  the  code  necessary  to  read  the  external  file  FT20F001  and 
manipulate  the  values  there  as  necessary  to  compute  the  desired  non-linear 
regression  covariance  estimates.  SAS  only  automatically  outputs  the 
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required  to  find  the  values  of  the  (asymptotic)  variance  and  covariance 
estimates. 

finally,  the  last  section,  beginning  with  the  statement 
PROC  PRINT  DATA  =  FINALOUT  LABEL; 

is  necessary  to  print  out  the  data  sets  in  which  all  the  estimates  have  been 
stored. 


5.  Output  Results 

(a)  For  the  linear  regression  procedure, the  output  results  consist  of 
the  following: 

-  Estimates:  Rx  =  ;  Ry  = 

-  Variance-covariance  of  Rx  and  Ry :  Var[Rx],  Var[Ry],  Cov[Rx  Ry]. 

-  Variance  of  In(P^)  =  (mean  squared  error) 

(b)  For  the  non-linear  regression  procedure  the  output  results  consist 
of  the  following: 

-  Estimates:  D0  =  ;  Rx  =  ;  Ry  = 

-  Variance  -  covariance  estimates  of  Dq,  Rx  and  Ry: 

VarCD03 ,  Var[Rx],  Var[Ry],  Cov[Rx,Ry],  Cov[Rx,D0],  Cov[Ry,r)0] 
-Variance  of  P^,  =  (mean  squared  error) 


Hardware:  IBM  4341  running  VM/SP  release  3  OASD  3370  (Disk  Drive) 
Operating  System:  VM/SP  CMS 

Programming  Language:  SAS79.6  with  SAS  ETS  and  SAS  GRAPH 


SECTION  IX 


COMPUTATIONAL  RESULTS 

In  this  section,  using  a  selected  set  of  13  matrices  for  four 

weapon/target  situations  at  different  impact  angles  ranging  from  15°  to  75°, 
at  estimates  of  Dg,  Rx,  and  Ry  are  given  using  the  non-linear  regression 
technique.  The  asymptotic  variance-covariance  matrix  of  Dg,  Rx,  and  Ry  is 
al so  given. 

For  reference  purpose,  the  result  of  the  linearized  regression  on  the 
two-parameter  Carleton  damage  function  (dg=l)  is  also  given.  This  includes 
the  estimates  of  Rx  and  Ry  and  the  variance-covariance  matrix  of  Rx  and  Ry. 
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1.  Estimates  of  Weapon  Parameters  1 


Weapon:  5EAKL  Target:  8T1.S  Impact 

Vel :  900  ft/sec 


Linear  Regression:  (Dq=1) 

-  Estimates:  Rx  =  233. 733  ft;  Ry  =  335.446  ft 

-  Variance-covariance  matrix  of  Rx  and  Ry: 


"x 

«y 


16.4307 

-4.8609 


31.0638 


-  Root  mean  square  error  on  In  Pj,:  /67295f?  =  2.5092 


Non-linear  Regression 

-  Estimates:  D0  *  .7551462;  Rx  =  72.9501  ft;  Ry  «  269. 2.381 Jt 

-  Variance  -  covariance  matrix  of  Og,  Rx> and  Ry: 


D0 


.000202 

.000250  1.3899 

-.000300  -1.7099  17.7441 


-  Root  mean  square  error  on  P^:  /.OOEOASSl  -  .04524 
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2.  Estimates  of  Weapon  Parameters  2 


Weapon  :  5EAKL 


Target :  8T1.5 

Vel :  900  ft /sec 


Impact  Angle:  30e 


Linear  Regression:  (Dg=l) 

-  Estimates:  Rx  =  214.054  ft;  Ry  =  343.765  ft 

-  Variance-covariance  matrix  of  Rx  and  Ry: 


Rx  14.1192 

R„  -5.3535 


26.5136 


-  Root  mean  square  error  on  In  P^:  /  5. 8167  -  2.4118 


Non-Linear  Regression 

-  Estimates:  00  =  .80309;  Rx  =  69. 3650  ft ;  Ry  =  292.5472  ft 

-  Variance  -  covariance  matrix  of  0g,  Rx, and  Ry: 


D0  .000227 


.000230  1.24858 


Ry  L  -.000330  -1.7550  20.8858 


-  Root  mean  square  error  on  P^:  00280587  *  0.05297 
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4.  Estimates  of  Weapon  Parameters  4 


Weapon  :  5EAKL  Target:  8T1.5  Impact  Angle: 

Vel :  900  ft /sec 


Linear  Regression:  (Dq=1) 

-  Estimates:  Rx  =  253.001  ft;  Ry  »  370.965  ft 

-  Variance  -  covariance  matrix  of  Rx  and  R  : 


Rx  9.53298 

Ry  -6.1764  23.1297 


-  Root  mean  square  error  on  In  Pk:  /3.369S  =  1.8357 


Non-Linear  Regression: 

-  Estimates:  D0  =  .85805  ;  1^  =  97.4700  ft;  Ry  =  317.1730  ft 

-  Variance  -  Covarlanve  matrix  of  DQ,  Rx, and  Ry: 

. 00027 

.00035  2.5862 

-.00038  -2.8052  25.8300 

-  Root  mean  square  error  on  P^:  /. 00418549  -  .0647 


60° 
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5.  Estimates  of  Weapon  Parameters  5 


Weapon  :  5EAKL 


Target:  8T1.5  Angle 

Vel :  900  ft /sec 


Linear  Regression:  (Dg=l) 

-  Estimates:  Rx  =  328.089  ft;  Ry  =*  362.266  ft; 

-  Variance  -  covariance  matrix  Rx  and  Ry: 


Kx  12.5726 

Ry  -4.0967  13.5213 


-  Root  mean  square  error  on  In  P^:  / .0105  =  .10247 


Non-Linear  Regression 

-  Estimates:  J)0  -  .92205;  ^  =  143.0987  ft;  Ry  =  366. 1414 jt 

-  Variance  -  covariance  matrix  of  DQ,  Rx,  and  Ry-. 

. 000306 

.000510  5.45596 

-.00040  -4.3990  30.21041 

-  Root  mean  square  error  on  P^:  / . 00743  =  . 086197 


D0 

RX 

Ry 
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6.  Estimates  of  Weapon  Parameters  6 


Weapon  :  5EAKL  Target:  US _  Impact  Angle: 

Vel :  900  ft /sec 


Linear  Regression:  (Dq=1) 

Ry  =  112.712  ft 
of  Rx  and  Ry: 

3.61066 

-.61024  4.28234 


-  Root  mean  square  error  on  In  P^: 

✓13.5915  =  3.6867 

Non-Linear  Regression 

-  Estimates:  Dn  =  .33256; 

Rx  -. 

12.9288  ft;  Ry  =  37.: 

-  Variance  -  covariance  matrix  of  1 

Dg,  Rx,  and  Ryi 

no 

.00092 

.0010 

.0100 

Ry 

-.0001 

-.0965  .0727 

-  Root  mean  square  see  error  on 

:  ✓. 00029504  =  .0172 

Estimates:  Rx  =  74. 5391  ft : 


Variance  -  covariance  matri) 


15° 
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7.  Estimates  of  Weapon  Parameters  7 


Weapon  :  5EAKL  Target:  16  Impact  Angl 

Vel :  900  ft/sec 


Linear  Regression:  (Dg=l) 

-  Estimates:  Rx  =  56.4879  ft;  Ry  =  118. 59  ft 

-  Variance  -  covariance  matrix  of  Rx  and  Ry: 


Rx  1.2879 

Ry  -.57349  3.55144 


-  Root  mean  square  error  on  In  P^:  710,7703  =  3.2818 


Non-linear  Regression 

-  Estimates:  Og  =  ,367^34;  Rx  =  11.4333  ft ;  Ry  =  41,9625  ft 

-  Variance  -  r.ovarifnca  matrix  of  Og,  Rx,and  Ry: 


DO 


.00285 
.  00239 
-.00293 


1.99746 

-.24438  23.5749 


-  Root  mean  sqji.  ?  error  on  :  /. 0105071  =  .10250 


8.  Estimates  of  Weapon  Parameters  8 


Weapon  :  5EAKL  Target:  16 _  Impact  Angle: 

Vel :  900  ft /sec 


Linear  Regression:  (Dq=1) 


-  Estimates: 

Rv  =  54.7351  ft 

;  Ru  =  123.717  ft 

-  Variance  - 

covariance  matrix  of  Rx  and  Ry: 

Rx 

.947699 

"y 

-.69615  3.14925 

_ 

-  Root  mean  square  error  on  In  :  /8. 62187  =  2.9363 


Non-Linear  Regression 

-  Estimates:  0g  =  .40749;  Rx  =  12. 1322  ft ; 

-  Variance  -  covariance  matrix  of  Dg,  Rx> and 


Ry  =  44.6990  ft 

V 


DO 


.00017 


Rx  .00014 

Ry  -.00017 


.11126 
-.  13663 


1.3250J 


-  Root  mean  square  error  on  P^.:  /. 0007 2634  -  .02695 


46° 
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9.  Estimates  of  Weapon  Parameters  9 


Weapon  :  5EAKL  Target:  16 _  Impact  Angle:  60 

Vel :  900  ft /sec 

Linear  Regression:  (Dg=l) 

-  Estimates:  Rx  =  71.0691  ft;  Ry  =  124.729  ft 

-  Variance  -  covariance  matrix  of  Rx  and  Ry: 

Rx  .964991 
Ry  -.  58275  2.93055 

-  Root  mean  square  error  on  In  P^:  78.  23459  =  2.8696 


Non-  Linear  Regression 

-  Estimates:  Dq  =  . 44089 ;  Rx  =  17.1892  ft ;  Ry  -  44.4199  ft 

-  Variance  -  covariance  matrix  of  Og,  Rx. and  Ry: 


D0 

Rx 

"y 


.00018 

.00019 

-.00017 

h- 


.1944 

-.1675  1.1318 


-  Root  mean  square  error  on  :  7700075406  -  .0275 
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10.  Estimates  of  Weapon  Parameters  10 


Weapon  :  5EAKL  Target:  310  Impact  Angle: 

Vel :  900  ft /sec 


linear  Regression:  (Dg=l) 

-  Estimates:  Rx  =  195.359  ft;  Ry  *  293.402  ft 

-  Variance  -  covariance  matrix  of  Rx  and  Ry 


Rx  11.0944 

Ry  -3. 4227 


22.8314 


-  Root  mean  square  error  on  In  P^:  /4. 25801  -  2.0688 


Non-Li  near  Regression 

-  Estimates:  00  =  .82377;  Rx  *  78.1014  ft;  Ry  *  296.6297  ft 

-  Variance  -  covariance  matrix  of  Dg,  Rx,  and  Ry: 

.000244 

.00022  1.6230 

-.00029  -  2.0473  22.2833 

-  Root  mean  square  error  on  :  / . 004 204 03  =  . 0648 


15° 
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11.  Estimates  of  Weapon  Parameters  11 


Weapon  :  5EAKL 


Target :  310  ^Pact  An91e: 

Vel :  900  ft /sec 


Linear  Regression:  (Dg=l) 

-  Estimates:  Rx  •  215.154  ft;  Ry  =  326.333ft 

-  Variance  -  covariance  matrix  of  Rx  and  Ry  ■ 


«y 


6.30579 

-4.2457 


17.9445 


-  Root  mean  square  error  on  In  :  /2. 36540  =  1«5380 


Non-Linear  Regression 

.  Estimates:  D0  =  .90123;  Rx  *  m.UMjt:  Ry  =  324^7570^ 

-  Variance  -  covariance  matrix  of  0g.  Rx, and  Ry: 

.00030 

.00038  4.0974 

-.00036  -3.6502  27.9007 

-  Root  mean  square  error  on  Pk :  / ."0'0M55'8.T  -  .09415 


Do 

Rx 

«y 


60° 
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Weapon  :  5EAKL 


12.  Estimates  of  Weapon  Parameters  12 


Target :  310 


Vel :  900  ft/sec 


Linear  Regression:  (0q=1) 

-  Estimates:  Rx  =  285.001  ft;  Ry  =  321.833  ft 

-  Variance  -  covariance  matrix  of  Rx  and  Ry: 


Impact  Angle:  75" 


Rx  9.25963 


Ry  -2.8505  9.96642 


Root  mean  square  error  on  In  Pk:  /j, 97828  *  1.40651 


Non- Linear  Regression 

-  Estimates:  D0  =  .94993:  Rx  =  163.7629  ft;  Ry  =  359. 1702  ft 

-  Variance  -  covariance  matrix  of  Oq,  Rx,and  Ry: 


00  .00035 


.00039  7.82936 


Ry  I -.00041  -5.7288  35.9268 


-  Root  mean  square  error  on  :  /. 01576  -  ,12554 
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13.  Estimates  of  Weapon  Parameters  13 


Example 

Consider  the  following  situation: 

Weapon  5EAKL,  Target  8T1.5,  Impact  angle  15°.  The  results  of  the  non-linear 
regression  analysis  are: 

D0  =  .75515;  Rx  =  72.9501  ft;  fty  -  269.2381  ft 

on  =  .01414;  <JD  =  1.  1784  ft;  aD  =  4.2124  ft 
U0  x  Ry 

The  two-standard  deviation  confidence  interval  on  the  weapon  parameters  will 
be 


0  ±  2  on  =  .75515  t  .02828 
u  u0 

R„  i  2  a.  =  72.9501  ±  2.3568  ft 
x  Rx 

\  ±  2  ap  *  269.2381  ±  8.4248  ft 

y 
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SECTION  X 


COMPUTATION  OF  ERROR  IN  Pk 


i 

i 

r 

\ 

t 

I 


Consider  again  the  three-parameter  Carleton  damage  function. 


pk  *  pk(°0,Rx*Ry) 

2  2 

=  Dg  exp  [ - Dq  {^2  +  ^)]  (53) 

Rx  Ry 

Given  estimates  of  Dg,  Rx,and  Ry,  say  Dg,  Rx,  and  Ry,  as  well  as  their 
vari ance-covariance  matrix.it  is  required  to  find  an  estimate  of  E[Pk]  and 
Var[Pk]  (which  is  a  measure  of  the  error  on  the  estimate  of  Pk),  at  a 
particular  target  location  (x,y). 

The  methodology  outlined  in  [6]  will  be  applied  here  using  the  Taylor's 
series  estimation  method. 


1.  Estimation  of  E[Pw] 

As  a  first  approximation  one  can  write 


E[Pk]  =  Pk(50,  Rx  ft) 


D0  exp  (-  00  + 


r 

x  y 


(54) 


i 


i 


U 


c 


> 


l ... 


L_ 
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L_ 


2.  Estimation  of  Var[Pwl 

Expanding  Pk  *  Pk(Dg,  Rx»  ty  as  a  Taylor's  series  about  the  point 
(00,  Rx,  Ry)  and  retaining  only  the  first  order  terms  result  in: 


Pk(n0,  Rx,  Ry)  S  Pk  (Oq.  Rx»  Ry) 

3  P 

+  {[}„  -  nj  (•?;— ) 

0  50*  V  *y 

3P 

+  <RX  *  V  CaT'i  n  u  r 

x  nn*  Rx’  y 
3  P 

*  "V  -  V  <*r>  v  »x,  «y 

Transposing  and  squaring  both  sides  of  (55)  yields 

cpk(o„.  rx  v  -  yo0, *yn2 


(55) 


7  3Pl  p 

-  <°o  -  6o>  Oo>  v 

9  3f\  9 

+  (RX  -  V  ^  g  R  R 
x  °0’  Rx*  Ky 

7  3Plr  7 

’'V-V2^  6  „  , 

*  u0*  V  y 

3p  ^  p 

+  2(nQ  -  DQ)(Rx  -  Rx)  [(ao-Kap^)] 

a  x  u0,  kx,  Ky 

3  P  3  P 

*  mo  ■  Do><*y  •  V  ll4)(5";\.  IS,.  5y 

3Pt  3Pl 

*  Z(R«  ’  '«u,» '  V  is,.  «y 


(55) 
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Taking  expectations  on  both  sides  of  (56),  one  obtains  as  a  first 
approximation: 


Var[P„]  -  V»r[D0]  (jA* 

o  (,0*  Kx*  y 


Or.  - 

+  Var[Rxl  (g^-)  _ 


x  Rq,  Rx,  Ry 


orv  ? 

*  Var[R  1  («£■)* 

y  y  R  ,  R 

J  0’  x’  x 

3P  3  P 

+  2  Cov[00.Rx]  [ [ [gR-)] 

o  *  R0*  Rx’  Ry 

3P.  3Pt 

+  2  cov[n0,Ry]  [(afi-Kinr)]-  D-  5 
o  y  00.  Rx,  Ry 

3  P  3  P 

+  2  Cov[Rx,RyJ  [('3R_)(3R~)] 

x  y  dq,  Rx,  Ry 


The  expressions  for  the  partial  derivatives  are  given  in  (46),  (47), and  (48) 
and  are  repeated  here: 


3Pk  ,2  ?  2  J> 

JrTs  -  no(~  *  J7)l  exP(-  no(^V  + 

0  R  R  R  R 


x  y 


x  y 


3  p  n4,  ?  ? 

x  Kx  \  Hy 


3Pk  ?  Dn  „2  v2 

3R“=  ?y  3  Gxpl-°nf3+3^ 
^  Ry  Rx  Ry 

x  y 


Substituting  (58),  (59),  (60),  and  (61)  in  (57)  yields 

Var[Pk]  =  (1-X)2  e'2X  Var[0Q]  +  4x4  e_2X  Var[Rx3 

X 

4-  4y4  e'2X  Var[Ry3 

Ry 

-2 

+  4x2  — r  (1-X)  e  ^  Cov(Oq,Rx] 

ft; 


B2 


+  4y2  ^  (1-X)  e'2X  Cov[00,Ry] 


¥ 

y 


3-2  e‘2X  Cov[Rx,Ry] 


Example 


+  8x  y  * 
oJ 

Rx  Ry 


(6: 


Given:  weapon  5EAKL,  target  8T1.5.  impact  angle  15°,  the  following  are 


sample  data  from  the  P^  matrix: 


' 1 

— 

97.2 

121.5 

37.9 

.8974 

.8764 

75.8 

.4263 

.3920 

Estimates  of  E[Pit]  using  Roger's  Ellipse  Program 


Estimates  of  E[Pi,1  using  Non-Linear  Regression  Procedure 


The  estimated  E[Pkrs  are  the  same  for  both  positive  and  negative  values  of 
the  x' s. 


Example 

For  the  following:  weapon  5EAKL,  target  8T1.5,  impact  angle  15°,  the 
matrix  provides  the  following  points: 


97.2  121.5 
37.9  .8974  .8764 
75.8  .4263  .3920 


To  find  Pk  at  x=50  ft  and  y=l00  ft,  an  interpolation  scheme  is  used.  At 
x-37.9  ft,  interpolation  yields  the  following  Pk  value  at  y= 1 00  ft. 


(100  -  97.2)  (.8964  -  .8974) 
pk  =  .8974  +  ^ - ri2i1r-  97-.T! - 


=  .8974  + 


(2.8)(-.021Q) 


=  .8950 


At  x=75.R  ft,  interpolation  yields  the  following  Pk  value  at  y=100  ft, 

r  ,  (100  -  97. 2) ( .3920  -  .4263) 

P*  =  *4263  +  - (121  .'5  -“772) 


=  .4263  ♦ 


[2.8) (-.0343) 


Interpolation  at  y=100  ft  yields  for  x=5Q  ft. 


Example 

For  the  following:  weapon  5EAKL,  target  8T1.5,  impact  angle  15°,  the 
matrix  provides  the  following  data  points: 


97.2 

121.6 

-37.9 

.2470 

.1770 

-75.8 

.0683 

.0555 

The  absolute  values  of  x  will  be  used  in  the  sequel.  To  find  Pk  at  x=50  ft 
and  y=100  ft,  an  interpolation  scheme  is  used.  At  x=37.9  ft,  interpolation 
yields  at  the  following  Pk  value  at  y=100  ft. 


At  x=75.8  ft,  interpolation  yields  the  following  value  at  y=100  ft 
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Interpolation  at  y=100  ft  yields  for  x  =  50  ft: 


P 


k  = 


.2389  + 


(50  -  37.9)  (.0668  -  .2389) 
(75.8  -  37.9) 


=  .2389  + 


(12.1)  (-.1721) 
(37^ 


.1840 


Note  the  extreme  differences  in  Pk  between  (50,100)  and  (-50,100). 


Example 

For  the  following:  weapon  5EAKL,  target  8T1.5,  impact  angle  15°,  Roger’s 
ellipse  program  yields 


D0  =  .79;  Rx  =  75.1046  ft;  Ry  =  2R8.0415  ft. 

The  EC P|c.3  value  for  x=50  ft.  and  y=100  ft  is 

E[Pk]  =  (.79)  exp  [(-.79  [ C 75^1046^  +  ^288.0415^  M 

=  (.79  (.640599)  =  .50607 

Example 

For  the  following:  weapon  5EAKL,  target  8T1.5,  impact  angle  15°,  the 
two-parameter  linearized  Carleton  damage  function  ( r)Q=  1 )  model  yields: 

Rx  =  233.733  ft,  and  Ry  =  335.446  ft. 

The  E[Pkl  value  for  x  =50  ft  and  y=100  ft  is 

2  2 

F[Pk]  =  exp  [-  (~  +  ^)] 
x  y 

t  r t  50  ,2  f  100 

=  exp  {-  [ (233.733)  +  (335.446)  )) 

=  exp(-  .13463)  =  .87404 
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3.  Example 

Consider  the  results  given  in  Section  IX. 1  for  the  following:  weapon 
5EAKL,  target  8T1.5,  impact  angle  15°. 

Here 

Dq  »  .7551462;  Rx  =  72.9501  ft;  Ry  =  269.2381  ft 

VarfOg]  =  .00020;  Var[Rx]  =  1.3899  ft2;  Var[Ry]  =  17.7441  ft2 

Cov[D0,Rx]  =  .00025;  Cov[D0,Ry]  =  -.00030;  Cov[Rx,Ry]  =  -1.7099 

It  is  required  to  determine  an  estimate  of  E[Pk]  with  a  two  standard 
deviation  confidence  interval  for  a  point  target  located  at  x=50  ft  and  y=100 
ft.  The  non-linear  regression  model  is  used  here. 

First  X  is  computed 

X  =  .7551462  U 7^950!  ]?  +  ( ^g^'Ml  ^ ^ 

=  .7551462  (.4697731874  +  .1379516711) 

=  .458921 

From  (54)  one  obtains 

W  =  00  e'X 

=  (.7551462)  e*- 458921 
=  .4772 
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From  (62)  Var[Pk1  is  computed 


Var[Pk]  =  (1  -  .458921) 2  e_2( ’458921)  (>n0020) 

♦  4(50)4  e-«-458M,>  (1.3OT9) 

(72.9501) 

♦  «U00)4  ( -7.5514621*  e-?(  .<58<I21) 

(269.2381) 

♦  4(50)2  (-7SS.14^^  (1  -  .458921)  e'2{ -45892l)  (.00025) 

( 72.9501 )  3 

+  4(100)2  (» 7551462 L2  {1  .  .458921)  e  ‘2<-458921)  (_.00030) 
(269.2381) 


+  8(50)2(100)2 


(.7551462) 

(7 2. 9501) 3 (2692381) 3 


-2(  .458921 )  (.K709g) 


Var[Pk)  =  .000023385  +  .000029942  +  .000002420  +  .000000794 
-  .000000140  -  .000005862 
=  .000.050,539 


The  standard  deviation  on  the  estimated  Pk  is 


'.000,050,539  = 


.0071 


A  two-standard  deviation  confidence  interval  on  Pk,  the  true  mean  value  of  Pk 

is 


Pk  =  ErPR]  ±  2op  =  .4772  i  .0142 

This  means  that  at  least  75  percent  of  the  future  values  of  Pk,  at  x=50  ft, 
y=100  ft,  should  lie  in  the  computed  interval. 
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SECTION  XI 


CONCLUSIONS  AND  RECOMMENDATIONS 

There  are  two  main  conclusions  that  can  be  reached  at  this  stage: 

1.  Roger's  ellipse  program  provides  a  good  approximation  to  the  estimate 
of  the  three  parameter  Carleton  damage  function.  These  estimates  compare 
favorab1y  with  the  estimates  obtaineo  from  a  non-linear  regression  procedure. 

2.  The  results  indicate  that  in  modeling  the  probability  of  kill  for  low 
impact  angles,  it  might  be  useful  to  consider  either  modifications  of  the 
Carleton  damage  function  or  other  functions.  One  particular  suggestion  is  to 
fit  two  Carleton  damage  functions,  one  for  positive  ranges  and  one  for 
negative  ranges. 

Some  of  the  problems  suggested  to  be  studied  in  future  work  are: 

1.  The  use  of  a  non-linear  regression  analysis  to  fit  a  blast  damage 
function  to  actual  data. 

2.  The  estimation  of  the  errors  in  when  fragmentation  anu  blast  are 
simultaneously  present. 

3.  The  development  of  a  methodology  to  identify  the  sources  of  errors 
and  to  measure  the  magnitudes  of  these  errors  which  are  inherently  present  in 
the  elements  of  the  matrix.  This  would  require  a  thorough  understanding 
and  analysis  of  the  methodology  used  in 

a.  the  fragmentation  testing  program  for  weapons; 

b.  the  target  vulnerability  program. 
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APPENDIX  A 

ROGER'S  ELLIPSE  PROGRAM 


C  ROGER'S  ELLIPSE  PROGRAM 
ASTH  -  0. 

DO  210  1-1,36 
RBSUM1  =■  0. 

R2SUM1  -  0. 

READ( 1 ,END=355)  (PK<J) ,J-1,NR) 

DO  200  J  -  1, NR 
RID  -  R1(J+1)-R1(J) 

R1S  -  R1(J+1)+R1(J) 

RBS  =  PK(J)*R1D 
RBSUM1  =  RBSUM1  +  RBS 
R2S  -  RBS  *  RlS 
R2SUM1  -  R2SUM1  +  R2S 

200  CONTINUE 

IF(R2SUMl.EQ.O.)  S2TH  =  1. 

IF(R2SUMl.EQ.O.)  GO  TO  201 
S2TH  -  1 . -RBSUM1 **2 /R2SUM1 

201  ASTH  =  ASTH  +  S2TH 
RTH  =  SQRT (R2SUM1 ) 

RT( I )  -  RTH 

STH  -  SQRT(S2TH) 

210  CONTINUE 
REWIND  1 

CONST  -  4.0/3,14159265 
AA  -  ASTH/ 36. 

DO  -  CONST* (l.-AA) 

NZ  -  36 

PI  =  3.14159265 
CA  -  0. 

Q  -  0. 

SA  =  1. 

DO  300  1  »  1,NZ 
AI  =  I 

TH  -  5. *(AI- . 5) 

TH1  =  TH/57.29578 
ST(I )  -  SIN(THl) 

ST2( I )  -  ST(I)**2 
CT(I)  -  COS(THl) 

CT2( I )  -  CT( I )**2 
300  CONTINUE 
320  D  =  0. 

DO  400  I  -  1,NZ 

TI  -  SQRT(CMAE/(PI*(SA*ST2(I)+CT2(I)/SA))) 
WI  -  TI**3*((SA*ST<I))**2-CT2(I))/SA**2 
RMT  «  RT(t)  -  Tl 
D  -  D  +  (RMT*WI) 
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400  CONTINUE 

IF(D.GT.O.)  GO  TO  309 
CA  =  SA 
SA  =  SA  +  1. 

GO  TO  320 

309  CONTINUE 

310  S  =  0. 

D  -  0. 

DO  500  I  =  1  ,NZ 

TL  =  SQRT(CMAE/(PI*(SA*ST2(I)+CT2(I)/SA))) 
W I  =  TI**3*((SA*ST(I))**2-CT2(I))/SA**2 
RMT  =  RT(I)  -  TI 
S  =  S  +  RMT**2 
D  =  D  +  (RMT*WI) 

500  CONTINUE 
S  =  2 . *S 

I F ( ABS (Q-SA) . LE . .001 )  GO  TO  5000 
I F ( D . GT . 0 . )  GO  TO  330 
CA  =  SA 
CO  TO  340 
330  C  =  SA 
340  Q  =  SA 

SA  -  (CA  +  Q)*-5 
GO  TO  310 
5000  SA  =  l./SA 

WRITE  (6,9)  SA 
RAD  =  SQRT(CMAE/PI) 

RSA  =  SQRT(SA) 

RX  =  RAD/RSA 
RY  =  RAD/RSA 
WKITE(6 , 350)  RAD 
WRITE(6 , 351 )  RSA 
WRITE(6 , 352 )  RX 
WKLT'ii(6,353)  RY 
360  FORMAT ( 1 HO , 8X , 'DO  =  '.F10.2) 

WRITE( 15 )  RX , RY , DO 
REWIND  l 
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APPENDIX  B 


SAS  COMPUTER  PROGRAM  FOR  ESTIMATING 
PARAMETER  VALUES  OF  PR 

* - * 

*PART  1 

THIS  SECTION  CONSISTS  OF: 

(1)  A  DATA  STEP  TO  READ  THE  INPUT  DATA,  GENERATE  THE  RANGE  AND 
DEFLECTION  VALUES  CORRESPONDING  TO  EACH  'PROBABILITY  OF  KILL' 
VALUE  (Z),  AND  CREATE  A  DATA  SET  WITH  THE  RANGE,  DEFLECTION, 

AND  Z  VALUES  EACH  STORED  IN  ONE  OBSERVATION.  I.E.  3  TOTAL 
OBSERVATIONS,  EACH  CONSISTING  OF  AN  IDENTIFICATION  (ID) 

VALUE  (X ,Y,OR  Z)  AND  1600  DATA  VALUES. 

(2)  A  PROC  STEP  TO  TRANSPOSE  THE  ABOVE  DATA  SET  INTO  A  DATA  SET 
WITH  1600  OBSERVATIONS,  EACH  CONTAINING  ONE  Z  VALUE  AND  THE 
CORRESPONDING  X  ANu  Y  VALUES. 

(3)  A  DATA  STEP  TO  ADD  THE  VALUES  LOG(Z)  (Z.GT.O),  X**?,  AND 
Y**2  TO  EACH  OBSERVATION.; 

* _ *•  • 

9 

DATA  TEMP , WPNDAT ; 

ARRAY  K  K1-K1600; 

INFILE  MYDATA; 

INPUT  IS  $  INT  JNUM  ; 

IF  ID-'Y'  THEN  DO; 

DO  1-1  TO  A; 

DO  J-I  TO  40; 

DO  L“1  TO  10; 

_Ij-L+40C*(I-l )+I0*( J-I ) ; 

K>(10*(I-l)+L)*INT; 

END;  END;  END; 

OUTPUT; 

END; 

IF  ID-'X'  THEN  DO; 

DO  1-1  TO  A; 

DO  J-l  TO  40; 

DO  L-I  TO  10; 

_I__-L+400*(I-1  )+10*(  J-L) ; 

K-(J-JNUM)*INT; 

END;  END;  END; 

OUTPUT; 

END; 

IF  ID-'Z'  THEN  DO; 

INPUT  K1-K1600; 

OUTPUT; 

END; 

RUN; 

PROC  TRANSPOSE  OUT-TPOSE; 

VAR  K1-K1600; 

ID  ID; 

RUN; 

DATA  LINPOSR; 

SET  TPOSE; 

XS-X*X; 

YS-Y*Y; 

IF  Z  GT  0  TIEN  LZ— LOG(Z); 

ELSE  LZ-.; 

RUN; 
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♦PART  2 

THIS  SECTION  CONSIST  OF! 

(1)  THE  "PftOC  REG'  PROCEDURE  STATEMENT  NECESSARY  TO  RUN  THE 
DESIRE  LINEAR  REGRESSION  ON  THE  FINAL  DATA  SET  CREATED  IN 
PART  1  ABOVE.  THE  'OUTEST-VCOVMAT'  OPTION  IN  THIS  STATEMENT 
CAUSES  THE  PROCEDURE  TO  OUTPUT  A  DATA  SET  CONTAINING  THE 
REGRESSION  PARAMETER  ESTIMATES. 

(2)  TWO  DATA  STEPS  WITH  STATEMENTS  NECESSARY  TO  CREATE  A  DATA 
SET  CONTAINING  THE  PROPERLY  LABELED  PARAMETER  ESTIMATES. 

(3)  A  THIRD  DATA  STEP  TO  PERFORM  THE  COMPUTATIONS  FOR  THE 
LINEAR  REGRESSION  ESTIMATES  OF  RX,  RY,  VAR(RX) ,  VAR(RY),  AND 
COV(RX,RY). 

(4)  A  FOURTH  DATA  STEP  TO  CREATE  A  1600  OBSERVATION  DATA  SET, 
EACH  OBSERVATION  CONTAINING  ONE  DATA  POINT  (X,  Y,  Z)  AND 
THE  RX  AND  RY  ESTIMATES  FOUND  ABOVE.; 


PROC  REG  OUTEST-VCOVMAT  COVOUT; 

MODEL  LZ-XS  YS/N01NT  NOPRINT; 

OUTPUT  OUT-WPNINFO  P-LZHAT; 

RUN; 

DATA  TX  TY  TP; 

SET  VCOVMAT; 

IF  _MODEL_- ' XS '  THEN  OUTPUT  TX; 

IF  _MODEL_- ' YS '  THEN  OUTPUT  TY; 

IF  _MODEL_- '  *  THEN  OUTPUT  TP; 

RUN; 

DATA  OUT; 

MERGE  TY  ( RENAME- (YS-VARBHAT ) ) 

TX  ( RENAME- (XS-VARAHAT  YS  -COVAB)) 

TP  (RENAME-(XS-AHAT  YS-BHAT)); 

RUN; 

DATA  FINALOUT(KEEP-RXHAT  RYHAT  SIGMASQ  VARRXHAT  VARRYHAT  COVRXRY); 

SET  OUT; 

RXHAT -SQRT ( I /AHAT ) ; 

RYHAT -SQRT ( 1 /BHAT ) ; 

S [GMAS0-_S IGMA_*  *2 ; 

VARRXHAT-VARAHAT/ (4*AHAT**3 ) ; 

COVRXRY-COVAB / (4 *(AHAT*BHAT )**l . 5 ) ; 

VARRYHAT-VARBHAT/(4*BHAT**3 ) ; 

CALL  SYMPUTCJX', RXHAT); 

CALL  SYMPUTCJY' .RYHAT); 

RUN; 

DATA  JPOSE(KEEP«X  Y  Z  RXHAT  RYHAT); 

IF  _N_-1  THEN  SET  FINALOUT; 

SET  LINPOSE; 

RUN; 
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*PART  3 

THIS  SECTION  CONSISTS  OF: 

THE  STATEMENTS  NECESSARY  TO  CALL  THE  PROC  NLIN  PROCEDURE, 
PRECEEDED  AND  FOLLOWED  BY  'PROC  PRINTTO'  STATEMENTS.  THE  'PROC 
PRINTTO'  STATEMENTS  CAUSE  THE  PROCEDURE  OUTPUT  FROM  PROC  NLIN  TO 
BE  WRITTEN  TO  AN  EXTERNAL  FILE  WHICH  IS  REFERENCED  LATER  IN  THE 
PROGRAM  WITH  AN  'INFILE'  STATEMENT.  THE  EXTERNAL  FILE  WILL  HAVE 
THE  NAME  'FT20F001'. 

NOTE:  THIS  IS  NECESSARY  SINCE  THE  PROC  NLIN  PROCEEDURE  OUTPUT 
INCLUDES  CORRELATION  BUT  NOT  COVERIANCE  ESTIMATES.  THE 
EXTERNAL  FILE  IS  READ  BACK  INTO  A  SAS  DATA  SET  (PART  4) 
AND  THE  COVARIANCE  ESTIMATES  ARE  COMPUTED  FROM  THE 
CORRELATION  ESTIMATES  AND  PARAMETER  STANDARD  ERRORS  WHICH 
ARE  PICKED  OUT  AND  LABELED.; 


PROC  PRINTTO  UNIT-20; 

RUN; 

PROC  NLIN  DATA-JPOSE; 

PARMS  RX-&JX 

DO-O-5  TO  1.3  BY  .05 
RY=& JY ; 

BOUNDS 

0<RX0900 

0<DO<«1.3 

0<RY<“900; 

TEMP-  -DO*(X**2/RX**2  +  Y**2/RY**2); 

MODEL  Z»DO*EXP (TEMP ) ; 

DER.RX=  (DO**2)*EXP(TEMP)*2*(X**2)/(RX**3) 
DER.RY-  (DO**2)*EXP(TEMP)*2*(Y**2 )/(RY**3) 
DER . DO-EXP (TEMP ) * ( 1+TEMP ) ; 

RUN; 

PROC  PRINTTO; 


* - - - * 

*PART  4 

THIS  SECTION  CONSISTS  OF: 

THE  DATA  STEPS  AND  PROC  STEPS  NECESSARY  TO  READ  IN  THE  EXTERNAL 
FILE  CREATED  IN  PART  3  ABOVE  AND  CREATE  A  FINAL  DATA  SET  WHICH 
CONTAINS  THE  NON-LINEAR  PARAMETER,  VARIANCE.  AND  COVARIANCE 
ESTIMATES-; 

* - - - - - * 

DATA  TWO; 

INFILE  FT20F001; 

INPUT  PID  $  ; 

IF  P ID= ' RX '  THEN  DO; 

INPUT  PI  ??  ; 

IF  Pl=.  THEN  DO; 

INPUT;  DELETE; 

END; 

ELSE  DO; 

INPUT  P2  P3;  OUTPUT; 

END; 

END; 

IF  PI D= ' RY '  THEN  DO; 

INPUT  P1-P3;  OUTPUT; 

END; 

IF  PID=' DO '  THEN  DO; 

INPUT  P1-P3;  OUTPUT; 

END; 

IF  PID=* RESIDUAL'  THEN  DO; 

INPUT  P1-P3;  OUTPUT; 

END; 

RUN; 

PROC  SORT  DATA=TWO; 

BY  PID; 

RUN ; 

DATA  BA(EEEP=RX  STDRX)  8B(KEEP=RY  STORY)  BC(KEEP“DO  STDDO) 

BD(KEEP=MSE)  BE ( KEE P=CORXD  CORXY)  BF(KEEP-CORYD) ; 

SET  TWO; 

BY  PID; 

IF  FIRST. PID  THEN  DO; 

IF  PID='RY'  THEN  DO; 

RY=P1 ; 

STDRY=P2; 

OUTPUT  BB; 

END; 

IF  PID=’RX’  THEN  DO; 

RX=Pl; 

STDRX=P2 ; 

OUTPUT  BA; 

END; 

IF  PID=' DO  1  THEN  DO; 

DO=Pl ; 

STDDQ=P2 ; 

OUTPUT  BC; 

l*v:n  . 

IF  PID=' RESIDUAL'  THEN  DO; 

MSE=P3; 

OUTPUT  BD; 

END; 


IF  LAST . PID  THEN  DO; 

IF  PID-'RX'  THEN  DO; 

C0RXD-P2 ; 

C0RXY-F3; 

OUTPUT  BE; 

END; 

IF  PID»'RY*  THEN  DO; 

CORYD=P2; 

OUTPUT  BF; 

END; 

END; 

RUN; 

DATA  ATWO; 

MERGE  BA  BB  BC  BD  BE  BF; 

RUN; 

DATA  AFINAL; 

SET  ATWO; 

VARRX»STDRX*STDRX; 

VARRY=STDRY*STDRY ; 

VARDO=STDDO*STDDO ; 

COVRXRY  =STDRX*STDRY*CORXY ; 

COVRXDO”STDRX*STDDO*CORXD ; 

COVRYDO“STDRY*STDDO*CORYD ; 

KEEP  VARRX  VARRY  VARDO  COVRXDO  COVRYDO  COVRXRY  MSE  RX  RY  DO; 
RUN; 


*PART  5 

THIS  SECTION  CONSISTS  OF: 

THE  ’PROC  PRINT'  PROCEEDURES  NECESSARY  TO  PRINT  OUT  THE  LINEAR 
AND  NON-LINEAR  ESTIMATES  IN  THE  PROPER  FORMAT.; 


PROC  PRINT  DATA=FINALOUT  LABEL; 

VAR  RXHAT  RYHAT  SIGMASQ  VARRX HAT  VARRYHAT  COVRXRY; 

TITLE1  LINEAR  REGRESSION  ESTIMATORS; 

TITLE2  MODEL...  -LOG(P)=X**2*(l/RX)**2+Y**2*(l/RY)**2; 

LABEL  RXHAT=RX  RYHAT =RY  SIGMAS  TERROR  VARIANCE 

VARRXHAT=VAR(RX)  VARRYHAT=VAR(RY)  COVRXRY«COV(RX,RY ) ; 

RUN; 

PROC  PRINT  DAT A=AF INAL  LABEL; 

VAR  DO  RX  RY  MSE  VARRX  VARRY  VARDO 
COVRXRY  COVRXDO  COVRYDO; 

TITLE 1  NON-LINEAR  REGRESSION  ESTIMATORS; 

Tl ILE2  MODEL . . . P-DO*EXP (-DO*(X**2/RX**2+Y**2 /RY**2 ) ; 

LABEL  VARRX* VAR (XX)  VARRY-VAR(RY)  VARDO-VAR(DO) 
MSE-ERROR_VARIANCE 

COVRXRY=COV(RX>RY)  COVRXDO=*COV (RX , DO )  COVRYDO =COV(RY , DO) ; 

RUN; 
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